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ABSTRACT 



\Q , 26 November 2008 

^p^l The effects of discreteness arising from the use of the iV-body method on the 

^ . accuracy of simulations of cosmological structure formation are not currently well un- 

' derstood. In the first part of this paper we discuss the essential question of how the 

C/3 , relevant parameters introduced by this discretisation should be extrapolated in con- 

vergence studies if the goal is to recover the Vlasov-Poisson limit. In the second part 
of the paper we study numerically, and with analytical methods we have developed 
, recently, the central issue of how finite particle density affects the precision of results 

^ ■ above the force smoothing scale. In particular we focus on the precision of results for 

\ the power spectrum at wavenumbers around and above the Nyquist wavenumber, in 

lO ■ simulations in which the force resolution is taken smaller than the initial interparticle 

spacing. Using simulations of identical theoretical initial conditions sampled on four 
different "pre-initial" configurations (three different Bravais lattices, and a glass) we 
• obtain a lower hound on the real discreteness error. With the guidance of our ana- 

' lytical results, which match extremely well this measured dispersion into the weakly 

OO , non-linear regime, and of further controlled tests for dependences on the relevant dis- 

■ creteness parameters, we establish with confidence that the measured dispersion is 

not contaminated either by finite box size effects or by subtle numerical effects. Our 
results show notably that, at wavenumbers below the Nyquist wavenumber, the dis- 
persion increases monotonically in time throughout the simulation, while the same 
^ \ is true above the Nyquist wavenumber once non-linearity sets in. For normalizations 

typical of cosmological simulations, we find lower bounds on errors at the Nyquist 
wavenumber of order of a percent, and larger above this scale. Our main conclusion is 
that the only way this error may be reduced below these levels at these physical scales, 
and indeed convergence to the physical limit firmly established, is by extrapolation, 
at fixed values of the other relevant parameters, to the regime in which the mean 
comoving interparticle distance becomes less than the force smoothing scale. 

Key words: Cosmology; N-body simulation; discreteness effects 



1 INTRODUCTION 

Dissipationless cosmological A^-body simulations aim to re- 
produce the clustering of dark matter in the universe, as- 
sumed to be in the form of a microscopic particle with ex- 
tremely weak non-gravitational interac tions (for reviews s ee 
e.g. IBertschinger. (,,19981: Bagia (2005i): lDolag et al.l l|2008t )'). 
In the absence of an analytical treatment of the strongly 
non-linear regime, these simulations have become increas- 
ingly central in extrapolating the predictions of the current 
"standard model" of cosmology to the corresponding scales. 



Many kinds of observations now probe directly or indirectly 
the distribution of dark matter at these scales, and will do 
so with greater precision in the coming years. The resultant 
need for precision in the theoretical results makes more nec- 
essary than ever a better understanding of these simulations. 
This paper concerns one potentially important source of er- 
ror which is currently still poorly understood: rather than 
evolving numerically the theoretical Vlasov-Poisson equa- 
tions describing the self-gravitating dark matter, simulations 
employ the A^-body method in which the matter is sampled 
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by "macro-particles". The errors introduced, i.e., the dif- 
ference between the results of the finite A'^ simulation and 
those in the theoretical model (which corresponds to an ap- 
propriate N ^ oo limit), are not understood. This is the 
discreteness problem in cosmological A^-body simulation. It 
is a problem which has received, given its potential impor- 
tance, a very modest amount of attention (see references be- 
low). Further the existing literature on the issue is marked 
by a considerable diversity in its conclusions, both quali- 
tatively and quantitatively. Given the ever more pressing 
need for robust control on the very considerable precision 
required of simulations — a goal of one percent precision i s 
now typically considere d (see e.g. lHuterer fc Takadal (|2005h : 
iMcDonald et al.l \200^ ) — it is an issue which deserves at- 
tention. 

In this paper we first give a brief review of the problem 
of discreteness in cosmological A^-body simulation. We both 
describe briefly previous work by other authors on the issue, 
as well as some recent work by ourselves and our collabora- 
tors, in which we have developed new analytical approaches 
to describe discreteness effects both in the initial conditions 
of simulations and in their early time evolution. As a starting 
point we attempt here to give a precise explicit formulation 
of the problem of discreteness. This distinguishes notably 
the problem from strictly numerical issues (e.g. about the 
agreement of codes using different summation techniques). 
We emphasize in particular the necessity to establish, before 
any discussion about the quantification of errors, precisely 
how numerical simulations should be extrapolated to ap- 
proach the desired theoretical limit. Our conclusion is simply 
that an appropriate such extrapolation is one which takes 
the interparticle spacing £ to zero, at fixed values of the other 
relevant discreteness parameters. Further such extrapolation 
should be done keeping fixed the initial conditions, which — 
given that simulations are performed in a finite periodic sys- 
tem — means using the same realization (and modes) of the 
initial theoretical power spectrum. 

After this introductory discussion we turn to a numeri- 
cal and analytical study of the issue. The goal of this study 
is to answer the more practical question of how small the 
interparticle distance £ must be to attain convergence of 
physically relevant quantities to a desired precision. In par- 
ticular we focus on the issue, which is at the centre of some 
controversy in the literature, of the accuracy of results, at 
scales around and below the initial interparticle distance, 
of simulations which use a force smoothing scale smaller 
than this latter scale. To attempt to resolve the question we 
use here a simple numerical method to isolate errors which 
manifestly must arise from discreteness, and which therefore 
give a lower bound on the discreteness errors. We focus here 
on the two point correlation properties of clustering, but 
the method we use can be extended to any other quantity 
(e.g. mass function, merger rates). The essential difference 
between our study and the few previous such attempts of 
this kind (see references below) is that we use, as mentioned 
above, also an analytical formalism which describes fully 
the measured discreteness effects at sufficiently early times. 
This allows us to "calibrate" the errors, in the sense that 
it allows us to establish without doubt that the measured 
quantities arise from physical discreteness effects, and not 
from the other possible sources of dispersion in our results 
(poor numerical convergence, or finite size effects). We then 



study the further evolution of these errors into the non-linear 
regime, allowing us to place with confidence a lower bound 
on the true discreteness error in this regime. 

More specifically our tests, and principal conclusions, 
are as follows. We consider A''-body simulations in an EdS 
universe of identical theoretical initial conditions, given by 
a random Gaussian realization of a power spectrum P{k) oc 
fc-^ In our central test the simulations differ only in the 
choice of "pre-initial" configuration, i.e., the point distri- 
bution chosen to represent the uniform universe prior to 
the application of the perturbations corresponding to the 
given theoretical initial conditions. The canonical choices in 
the literature are "grid" (a simple cubic lattice), or "glass" 
l|Whitj|l993l '). Here we consider a wider class of such config- 
urations, employing also two different Bravais lattice config- 
urations (body-centred cubic and face-centred cubic). The 
reason for our choice of these configurations is that they al- 
low us to apply in a very powerful way our analytical treat- 
ment. This formalism gives a very accurate description for 
the early time evolution of simulations starting, in principle, 
from any perturbed Bravais lattice. The differences in the 
evolved power spectra, starting from the same realization 
of the theoretical model discretised on these different distri- 
butions, measured in our simulations are in extremely good 
agreement with these analytical predictions for all wavenum- 
bers at early times, and progressively deteriorates, as antic- 
ipated, as we go into the strongly non-linear regime. In this 
latter regime we observe that these differences show similar 
dependences on the discreteness parameters as in the regime 
where we can fit them analytically. Notably at a given phys- 
ical scale, they decrease as £ decreases, and increase mono- 
tonically as a function of time at fixed £. 

These tests give us a robust non-trivial lower bound 
on the size of systematic discreteness errors. For the power 
spectrum, which is the quantity we focus on, these lower 
bounds are of order a few percent for wavenumbers compa- 
rable to and larger than the Nyquist frequency for a start- 
ing red-shift equal to 2^, and then decrease monotonically 
at smaller wavenumbers. While the precise bounds for any 
given cosmological model (and choice of other relevant sim- 
ulation parameters) will differ, they will be of this order (or 
larger, as the bounds monotonically increase with the start- 
ing redshift). Our results allow us then to draw conclusions 
about the question of how far £ must be extrapolated to 
attain errors smaller than of this order. Specifically we con- 
clude that the common practice of using results considerably 
below the scale £ (or tt/^ in reciprocal space) may be justi- 
fied, but only with a discreteness error bar which is, for the 
power spectrum, and for normalisations typical of cosmolog- 
ical simulations, of order of these lower bounds, i.e., several 
percent. Precision greater than this for a given wavenumber 
k, e.g., down to below the one percent level for the power 
spectrum now often cited as necessary, can be achieved only 
by using particle densities such that k£ < 1. Indeed results 
of simulations do not converge to the continuum limit until 
this parameter range is reached, and thus one cannot have 
real confidence in results without performing such an ex- 
trapolation. While this conclusion has been argued for in 
several studies by some authors (see references and discus- 
sion below), it is a much more stringent requirement than 
that assumed in much of the literature, and formulates a 
considerable challenge to simulation. 
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2 THE PROBLEM OF DISCRETENESS 

The problem of discreteness in cosmological simulation 
arises from the fact that the numerical simulations are not 
a direct discretisation of the equations of motion of the 
theoretical model. The latter is (usually) assumed to be 
described, on the physical scales of relevance, by Vlasov- 
Poisson (VP) equations (or "collisionless Boltzmann equa- 
tions" ) which give the evolution of the (smooth) phase space 
mass density. A'^-body simulations, on the other hand, are 
numerical integrations of the equations of motion of N self- 
gravitating particles, i.e., 

Gm . 

n3 



a 

2-x, 

a 



TV,(|x.-x,|) (1) 



where dots denote derivatives with respect to time, m is the 
mass of the particles, Xi is the (comoving) position of the i-th 
particle, and We(|xi— Xjj) isafunction which regularizes the 
divergence in the gravitational force at |xi — Xj| =0 below 
a characteristic scale e. These unphysical "macro-particles" 
are artefacts of the A''-body simulation technique, with a 
mass many orders of magnitude (typically ~ 10^") larger 
than those of the theoretical dark matter particles. 

As the VP equations may, in principle, be obtained as 
an appropriate N ^ oo limit of the particle system, the 
problem of discreteness is in practice that of determining the 
discrepancy between the solution of the A'^-body equations 
for some finite A'^ and their solution for a much much larger 
A'', representative of the VP limit. It is therefore evidently 
essential to specify precisely how to extrapolate cosmological 
A'^-body simulations to this limit. This is the point we first 
discuss. 



2.1 Discreteness parameters 

In the case of the A^-body method employed to solve the 
cosmological problem, the unphysical parameters character- 
ising the numerical solution can be clearly divided into two. 
Firstly there are those required, in addition to the parame- 
ters of the input theoretical model, to characterise the equa- 
tions lU and their initial conditions. Secondly there are the 
parameters introduced to then solve these well posed equa- 
tions numerically (e.g. time step, parameters controlling the 
precision of the calculation of the force). It is only the for- 
mer, which we will refer to as the discreteness parameters 
and denote by {Da}, which are the subject of study here. 
The latter set of parameters, which we will refer to as the 
numerical parameters of a simulation, control the accuracy 
with which the set of equations ((T}, with well defined ini- 
tial conditions, are solved. They therefore have no relevance 
to the problem of discreteness which we are focussing on: 
we wish to understand the relation between the results of 
a "perfect" A''-body simulation, i.e., an arbitrarily precise 
numerical solution of the equations ([1} from well specified 
initial conditions, and the evolution of the theoretical model 
from its corresponding initial condition^ 



^ The sensitivity of results to this second set is, of course, essen- 
tial to understand in order to characterise the precision of results 
of simulations (e.g. using different codes), and indeed considerable 
effo rt to improve con t rol ha s been made i n the l ast few years (see 
e.g. iHeitmann et al.l ||2007| '): iLukic et all ll2007t) '). We note that 



The set of discreteness parameters {Da } we consider is 
the following: 

• 1. The mass of the macro-particles (referred to as mass 
resolution in the literature), or equivalently (since the mean 
mass density is specified) their mean (comoving) number 

X /3 

density no . We will parametrize this by ^ = rig , which we 
refer to as the mean interparticle spacing. 

• 2. The smoothing parameter e characterising the regu- 
larisation of the force (known as the force resolution in the 
Uteratur^EI) • 

• 3. The pre-initial configuration, i.e., choice of grid, glass, 
or other distribution. We will denote this discrete variable 
prelC. 

• 4. The initial red-shift, Zi. 

Some remarks on this list are appropriate: 

• Discreteness effects depend on the number of particles 
A'' used in a simulation only through the density of these 
particles. Change in results when A'' varies, at fixed £, is not 
a discreteness effect: to pass from the VP equations to the 
TV-body equations ([1]), and their initial conditions, we do not 
need to introduce the side L of the cubic box (which, with 
periodic replicas, is canonically used to approximate the in- 
finite universe) as both sets of equations are well defined in 
infinite space. The box size L thus belongs to the second set 
of parameters, as it is introduced to solve the Eqs. ^ in 
a (finite) numerical simulation. The dependences on it, i.e. 
on the variation of A'' at fixed particle density, are finite-size 
effects. We do not study these effects here, and will always 
work at fixed L in our numerical study below. For studies 
of them see e.g. IPeiil Il997h : ISirkol (|2005l ): iBagla fc Prasad! 
(|2007l ): lBagla et al.l (2008) 

• The smoothing parameter e, on the other hand, cannot, 
in modern cosmological simulation, be considered as belong- 
ing to the numerical parameters: it is not, in this context, a 
parameter introduced to facilitate the numerical solutior(f|. 
Rather, as we will discuss further below, it is used with the 
aim of reducing efi^ects of two-body coUisionality, i.e., to try 
to make the TV-body solution approach better the theoretical 
collisionless behaviour corresponding to the VP equations. 

• That the initial red-shift Zi is a discreteness parameter 
in the sense we have d efined above has bee n shown explicitly 
inlJovce et al.l (|2005h : lMarcos et al.l l|2006l ): |jovce fc Marco^ 
(l2007al 'l fand summarised also briefly in Sect. l2.4l belowV Put 
simply, the treatment of the evolution of Eqs. ([TJ in this 
work shows, analytically, that the initial conditions for the 
TV-body system derived for a given input power spectrum 
(PS) at a redshift zi (using the canonical method based on 
the Zeldovich approximation), do not evolve exactly under 
Eqs. IT} to those set up from the same PS at a different red- 
shift Z2. This is true in the limit of arbitrarily small initial 

the distinction we make here between the two kinds of parame- 
ters is not usually m ade in the literatur e on the "convergenc e" of 
simulations (see e.e. lPower et al.l ||2003^ : iLukic et al.l fcOOTl )'). 
^ It is often referred to simply as the "spatial resolution" . We will 
not use this nomenclature here as the central issue we discuss is 
whether such an identification of the force smoothing scale with 
that of the spatial resolution is valid. 

^ Indeed the equations JTJ may be solved numerically without 
any such smoothing (and often are in other contexts e.g. galactic 
dynamics). 
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relative displacements to the lattice where non-linear (fluid) 
corrections to the Zeldovich approximation can be neglected. 

This list of discreteness parameters is a minimal one, appro- 
priate for, say, a standard P^M type code. Even in this case 
it could be elaborated to be more precise. For example, the 
regularisation involves the choice of a function which is not 
always the same, and e can vary in timsQ. Different choices 
of the sampled modes may also be made in setting up ini- 
tial conditions. The list is adequate also for a simple PM 
code, but evidently it would need to be expanded to describe 
adaptive codes in which the particle number changes in time 
and space according to some criteria. We will not consider 
such complexities here, apart from a few further comments 
on this point in our conclusions: it is sufflciently ambitious 
to hope, at least as a first step, to control fully the effects of 
discreteness for these simpler cases. 

2.2 Convergence to Vlasov-Poisson limit 

Let us now denote by Q{r, z; {Da}) the measured value of 
any physically relevant quantity in an A'^-body simulation, at 
red-shift z, with values of the discreteness parameters Da, 
e.g., a two point correlation function or a PS (where the 
variable r is given the appropriate interpretation, and could 
equally well represent a set of vector separations for a higher 
order statistic). 

The discreteness problem can be schematically repre- 
sented then as that of determining an estimate of the differ- 
ence 

AQ(r,z;{Da}) = Q{r,z;{Da}) - Qvpir,z) (2) 

where Q{r, z; {Da}) is the result of a "perfect" A^-body sim- 
ulation, and Qvp{r,z) is the result of the same quantity in 
the VP equations evolved from the same initial conditions. 

By construction Qvpi^jZ) is, in general, unknown. In- 
deed it is because we cannot determine it analytically that 
we turn to A'^-body simulation. To estimate AQ(r, z; {Da}) 
the best one can do is thus to study, numerically, the conver- 
gence of Q{r, z; {Da}) towards some fixed value as the {Da} 
are appropriately extrapolated. If the goal is to approach 
as closely as possible the evolution of the VP equations, 
one should evidently extrapolate the relevant parameters in 
a way which indeed gives convergence to this limit of the 
A'^-body system. While it is evident that the interparticle 
distance £ should be decreased, how the other parameters 
should be varied (or not) is not. Indeed, as we will discuss 
further in our conclusions, most of the few convergence stud- 
ies of cosmological simulations in the literature do not adopt 
an extrapolation which converges directly to the VP limillj. 

There is in fact no rigorous treatment in the literature 
on cosmological A'^-body simulations, or more broadly in 
the cosmology literature, establishing the existenc e of the 
VP limit: deriv a tions of the VP equations (see e.g. IPeeblesI 
1 19801 ): ISaslawl (|l989i ')') are limited to showing that these 

* We have implicitly assumed it to be fixed in comoving length 
units, which is usually the case, although many other variants can 

be found in the literature. 

^ A n otable exception is the study reported in ISplinter et al.l 
lll998l'). which w i ll we d iscuss below, as well as a recent paper 
bv lRomeo et al.l ||2008i). 



equations may be obtained by a truncation to the leading 
term of a BBGKY hierarchy of equations, but do not rigor- 
ously establish the conditions under which the required trun- 
cation may be mad^f]. Formal proofs establishing the validity 
of the Vlasov mean field approximation for long-range inter- 
acting systems can, however, be found in t he math e matic al 
physics literature (for a discu ssion see e.g. ISpohnI lIlQQll )). 
Notably iBraun fc HeppI (|l977 ') have proved that in a finite 
system of particles interacting through pair forces, reg- 
ularized so that the potential is bounded below at r = 0, the 
Vlasov limit corresponds to A'^ — > oo. In taking this limit 
the volume, mass and time of evolution are kept flxecQ. We 
will assume, without rigorous proof, the evident extension 
of this result to the inflnite volume case of cosmological sim- 
ulations: we take the VP limit as ^ ^ (i.e. particle number 
in any finite volume goes to infinity) at fixed mass density, 
followed by £ ^ The convergence at fixed temporal dura- 
tion corresponds to keeping also the initial red-shift Zi fixed, 
and as the limit should clearly not depend on preIC (the 
choice of pre- initial configuration), we also keep this fixed. 

In summary, applied to an A'^-body simulation, this tells 
us that an appropriate extrapolation is given by decreasing £ 
(i.e. increasing the particle density) while keeping the other 
discreteness parameters {Da} fixed. The limit is taken at 
fixed e 7^ 0, which means that the spatial resolution (for 
unsmoothed gravity) is limited to above this scale. In other 
words this extrapolation converges to a smoothed version of 
the VP equations, which then (we assume) would converge 
to VP as £ ^ 0. This extrapolation is not necessarily unique 
— convergence may in principle be obtained while allowing 
Zi and/or e to vary in various manners as a function of ^ — 
but it is certainly simple. The use of any alternative (if, we 
emphasize, the goal is to obtain direct convergence to the VP 
limit) should, however, be carefully considered to establish 
(at least as rigorously as here) that it gives convergence to 
the VP limit. 

It is important to note the specific order of the limits in 
£ and e. Beyond the necessity to introduce a regularisation in 
rigorous proofs of the VP limit mentioned above, the reason 
for this can be understood easily on physical grounds: the 
scale e 7^ provides a characteristic scale which is clearly 
necessary to give physical significance to the limit £ — > 0. In- 
deed taking e = and initial conditions specified by a pure 
power law input PS, £ is the sole characteristic scale of the 
discrete system (in the limit L — > oo), and defines itself the 
unit of length. Varying £ gives, up to a trivial rescaling, a 
system with exactly the same dynamics, which is manifestly 
not that of the VP limit (as it includes explicitly non-mean 
field effects such as two body coUisionality) . When, on the 
other hand, the same system is treated, but now with e 7^ 0, 



^ For an alternative derivation of the VP equations using a 
coarse-g raining of the micr oscopic equations for the particle sys- 
tem, see lBuchert &: Domingucz, (,,2005; ). 

^ Formally the coupling in the interaction (i.e. Gm^ for gravity, 
where m is the particle mass) scales in p r oport ion to 
* As noted, the proof of iBrauri fc Hepd 1 I1977I) is for finite non- 
zero smoothing, and the existence of the exact VP limit e — > has 
not in fact been proven. We neglect this mathematical subtlety 
here, which one would expect to be relevant, at most, to the 
asymptotically long time behaviour of the system (which does 
not interest us in this context) 
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the limit I Q has a non-trivial meaning if it is taken at 
constant e. Indeed non mean-field effects such as two body 
scattering, and other ones we will describe below, are explic- 
itly no longer present in the dynamics in this limit. 

Finally we make one important remark about initial 
conditions. Evidently to establish convergence to the evolu- 
tion of a given continuum model, one should keep the initial 
conditions in this limit fixed. While we have stated explicitly 
that, in extrapolating, the initial redshift Zi should be held 
fixed, this does not prescribe unambiguously the initial con- 
ditions at any finite I: given that simulations are performed 
in a finite volume, the number of modes in any interval of 
wavenumber is finite, and thus different realizations of the 
same initial conditions introduce intrinsic statistical fluctu- 
ations in the initial conditions compared to the average the- 
oretical behaviour. For convergence studies of discreteness 
effects specifically it is simplest to keep also the realization 
fixed, although of course such effects can in principle be av- 
eraged out by a sufficiently large number of realization^. 

2.3 How far must one extrapolate? 

While the above discussion simply tells us how to extrapo- 
late towards the VP limit, the practical form of the discrete- 
ness problem is more detailed: How small do we need I to 
he, given certain fixed values of other parameters in the set 
{Da}, to attain a desired precision AQ(r, z; {Da}) on the 
theoretical quantity Qvp{r,z) ? 

2.3.1 Common wisdom 

Current practice in interpreting the results of cosmological 
simulations appears to repose on an approximate answer to 
this question, which we now attempt t o summarize (see e.g. 
ISmith et all (|2003l '): |Power et all HqoI)). It is supposed that 
there are essentially two ways in which discreteness can play 
a significant role in making an A'^-body simulation deviate 
from the desired VP evolution: 

• Dl. Through the hmits placed on the accurate repre- 
sentation of the initial conditions. Indeed, to avoid alias- 
ing effects, only modes of the input theoretical PS up to 
the Nyquist frequency fcjv — tt/I should be sampled. Unless 
kc <C fcjv this means that there is "missing power". In sim- 
ulations of CDM type models, notably, this is always the 
case. Further, for any initial PS, there is always additional 
power in the initial conditions, predominantly at k > fcjv, 
generated purely by the discretisatiorF°l. 

^ Further, when £ > e, one needs to specify whether power should 
be added (if present in the theoretical model) in the larger range 
of wavenumbers which can be sampled in the initial conditions 
as £ decreases. While this is not the source of ambiguity in our 
extrapolation (as power should evidently not be added in the 
range that £ < e), for studies in the range £ > e, such as that 
we will report below, it is relevant. We will use the prescription 
that the realization of the input displacement field is kept fixed, 
as this allows us most clearly to identify effects of discreteness. 

Analytical expressions f or the full initial conditions are given 
in Ijovce fc MarcosI l l2007bl ). On a lattice, at linear order in an 
expansion in the amplitude of the input spectrum, the discrete 
power is non-zero only for k > fejv- lu a glass there is also a 
contribution, oc at small k, for k < fcjv. 



• D2. Through two body collisions in the course of the 
dynamical evolution which cause deviations from the desired 
mean field behaviour of the VP limit. 

The flrst point is believed not to place, in practice, an 
important limitation on the accuracy of simulations once 
they are evolved. The reason is that gravitational clustering, 
from CDM cosmological initial conditions, is understood to 
develop essentially by the transfer of power from large to 
small scales, non-linear structures being formed by the evo- 
lution of fluctuations at initially larger scales. The spatial 
resolution thus improves rapidly as time goes on, essentially 
following the forming non-linear structures which depend 
only on the presence of the initial fluctuations which seeded 
thenO. Small residual effects are envisaged, arising from 
the "spurious" power generated by the sampling on a spe- 
ciflc pre- initial configuration (grid or glass), but they are 
usually assumed to be negligible and of no practical impor- 
tancfO. An exception is in the case of hot (or warm) dark 
matter spectra. In this case one may have kc <^ kN so that 
all the initial power is well represented, but the small scale 
power generated by discreteness can evolve to form struc- 
tures which may not be "wiped out" sufficiently rapidly by 
the structures forming at larger scales. Recently interest in 
this case has been regenerated in the context of simulation 
of "warm dark matter" mod els, and it has been shown ex- 
plicitly in numerical studies l|Gotz fc Sommer-Larsenll2003l : 
IWang fc Whit3l2007l ). using different preIC (grid or glass) 
that such effects may be important, leading to gross dis- 
creteness effects in such simulations. 

The effects of two body coUisionality (D2) are under- 
stood to be taken care of by the smoothing e. Indeed it 
is explicitly for this reason that such a smoothing is intro- 
duced, its value being chosen ideally large enough to sup- 
press the related effects, but small enough so that too much 
spatial r esolution is not lost. Since , according to simple es- 
timates l|Binnev fc Tremainelll993 ). one expects such effects 
to be largest in regions of highest density, e is chosen just 
large enough to suppress them, over the relevant cosmolog- 
ical time scales, in such regions. 

In summary, these physical arguments may be formu- 
lated as qualitative answers to the question posed above, as 
follows. For typical quantities measured in simulations (e.g. 
two point correlation function, PS, halo masses and pro- 
files), the errors AQ(r, z\ {Da}) due to discreteness, for any 
r a little larger than e are negligible, i.e., so small as to be 
of no practical interest (compared to attainable numerical 
errors, notably), if: 

• Al. I is sufficiently small so that, at red-shift z, the 
fiuctuations at scale r may be formed by the collapse of 
fiuctuations initially at scales k < kN. 

• A2. I is sufficiently small so that the coUisional relax- 
ation time scale in the densest resolved regions (i.e. the high- 

A numerical stud y which nicely illustrating this may be f ound 
inlLittle et al.l lll99ll'). S ee also lBagla fc Padmanabhaiil lll997tl and 
iBagla fc PrasadI ||2008|). 

some works fe.E. ISmith et al.l 1I2OO3I') ') attempt to correct for 
the associated effects by subtracting this power which can be 
measured in the initial conditions. This procedure assumes that 
this spurious power does not evolve, a n assumption wh i ch we h ave 
shown analytically to be incorrect in I Joyce fc MarcosI (l2007bl) . 
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est density in a region of radius ^ e) is large compared to 
the age of the universe. 

Both of the answers can be converted, by making 
use of known phenomenological models describing the re- 
sults of simulation s (not ably halo models, or the model of 
iPeacock fc Dodd3 (|l996l )). into approximate criterion for 
the necessary £ (i.e. particle density) expressed in terms of 
the parameters of the theoretic al model, and the sca les r 
and £, and of the red- shift z. In iHamana et al.l (|2002l ) one 
can find, for exa mple, approximate c r iteria derived using 
hal o models, while [Knebe et al.l l|2000l ). iPower et all (|2003i ) 
and| piemand et al.l (|2004| ) present extensive numerical stud- 
le A series of other articles focus specifically on the ef- 
fects of two body relaxation in placing limits on the accu- 
racy of de nsity profiles in ha l os, us ir ig mostly numer i cal ap - 
proache s (iBirinev fc Knebd (120011 ). jPiemand et all ||2004) . 
lEl-Zand (|2006l 'l). A recent paper bv lBagla fc Prasad! ([20081 ') 
concludes, on the basis of some simple numerical tests on dif- 
ferent theoretical initial conditions, that discreteness effects 
may be neglected once the non-linearity scale has evolved to 
be larger than the mean interparticle separation. 



2.3.2 Dissenting views 

While these answers may be correct, they are certainly not 
in any way rigorous. The essential problem is that they as- 
sume that the physical effects of discreteness are known, or, 
at least, that those which play any significant role in simula- 
tions are known. While the latter may a posteriori prove to 
be true, the former certainly is not. Indeed understanding 
of the role of discreteness in the highly non-linear evolution 
of these systems is extremely limited. 

One of the surprising aspects, at least at first sight, of 
the standard criteria just discussed is that they allow the 
resolution scale of a simulation (at 2; = 0) to be very much 
smaller than the scale £. Indeed in practice the spatial res- 
olution is usually taken to be fixed by e, with e <^ £ 0. If 
one considers that this smoothing is introduced to make the 
"macro-particles" behave like fluid elements, moving under 
the effect of the mean field, it would appear to be neces- 
sary to have, at least, e ~ I. This point has been forcefully 
argued b y Melott et a l. in a series of pape r s during the 
nineties ([Melottl 1 199(3: iKuhlman et al.l 1 19961 : iMelott et al] 
ll997l:ISplinter et al.lll998| ) . and r estated in a recent com ment 
(lMelottil2007l ) InlMelottI (|l99(]| ). lKuhlman et all (|l996l ) and 
iMelott et all(|l997f )" specific non-Vlasov effects are explicitly 
shown to be present in numerical experiments, and of much 
greater importance once the regime e < ^ is attained. One 
of the few controlled studies of the issue of discreteness in 



In these latter papers the question of discreteness is not sepa- 
rated from the question of numerical convergence of the Af-body 
equation. Thus particle density and the smoothing e are consid- 
ered on the same footing as choice of time step, and parameters 
for force precision etc.. The sets of simulations studied do not in 
fact define a convergence study to the VP limit, as we have dis- 
cussed above: power is added as the particle density is increased, 
and the initial red-shift also changes. We will return to this point 
in our conclusions. 

In the "Millenium" simulation ISpringel et al. for ex- 

ample, £ fa <?/50. 



the liter ature in a spirit resem bling that advocated above is 
given in lSplinter etahl (|l998l ). The paper focuses on the dif- 
ference between results of simulation using PM codes and 
P^M codes at different resolutions. Its conclusion is that 
results of the latter codes in the regime e < I, do not agree 
well, most notably for phase sensitive statistics, with those 
obtained from higher resolution PM codes (for which e = I). 
These results, which place a question mark over the reliabil- 
ity of results below the scale I, have been largely ignored and 
addresse d only very incomp le tely in subsequ e nt wo rks (see, 
notably, iKnebe et ahl (|2000l ): iHamana etahl (|2002l )') which 
support, broadly, the "common wisdom" which we have out- 
lined above 

The common wisdom has al so been questioned b y 
several other groups of authors JSuisalu fc Saar (Il995h: 
iBaertschiger et al.1 (|2002l '): IXiao et al.1 (|2006l '): iRomeo et al.1 
( 2008al )). all placing in question (like Melott et al.) the use 
of a smoothing e < ^ on the basis of numerical results. In 
particular we note that the role played by interactions of par- 
ticles with their nearest neighbours — which give physical 
effects clearly not representative of the mean field Vlasov- 
Poisson limit — in the evolution of clustering at early times 
in simul ations has been highlight ed in cosmological simula- 
tions in IBaertschiger et al.1 (|2002l ). and i n a simplified class 
of grav it ationa l N-body simulations in IBaertschiger et ahl 
(|2007al lbl. l200i '). In a very recent studv lRomeo et al.1 (|2008ah 
conclude, on the basis of a study using wavelet techniques 
to analyse a set of ACDM simulations, also that results be- 
low the scale I are un reliable. We no te also the discussion of 
discreteness effects in iBinnevI (|2004l ). which illustrates with 
a study of a one- dimensional sheet model that discreteness 
may induce effects prior to virialisation (and distinct from 
two body effects) by artificially bounding above the growth 
of the phase space density. 



2.4 Analytical results 

In recent work b y ourselves and o ur collaborators 
(iJovce et ai] l2005l: [Marcos et ah] I2OO6I : iJovce fc Marcos! 
l2007al : iMarcosI '2008). we have used a perturbative treat- 
ment of cosmological A'^-body simulations to treat discrete- 
ness effects analytically. While the method is limited by its 
range of application (to sufficiently early times) it has the 
advantage of providing an exact quantification of these ef- 
fects in that range, as well as an understanding of the phys- 
ical mechanism at play. In this section we will briefly review 



We say "incompletely" because no other published work has, 
to our knowledg e, reported s i milar precise tests measuring the 
same quantities. iKnebe et al.l 1 I2OOCII) ascrib e the differenc e s seen 
in the two point correlation properties by ISplinter" et al.l (Il998l ) 
to "erroneous evolution in high resolution runs" , b ut without any 
roof (their own numerical tests, unlike those of iSplinter et al.l 
19981 ). are not tests for discreteness effects but for the coher- 
ence of result s prod uced by different codes). The analysis of 
iHamana et al.l 1I2OO2I) . which explicitly calculates how resolution 
improves with time as foreseen by the "common wisdom" de- 
scribed above, suffers from the weakness, underlined by the au- 
thors themselves, that it is based on the use of a halo model 
description of non-linear clustering, itself drawn from numerical 
simulations. 
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this formalism, which we will employ in the next section in 
the analysis of our numerical results. 

The treatment can be understood as a generalization 
to discrete distributions of the standard linearization of the 
equations of a se lf-Rravitating fl uid, in the Lagrangian for- 
malism (see e.g. iBuchertI l|l992l )). At linear order we thus 
refer to it as "particle linear theory" (PLT). The full de- 
tails can be found in these publications, and we will limit 
ourselves to a short summary of t he essential idea, and the 
salient results. We note that while [jovce fc Marcos! (|2007al ) 
presents the details of the use of PLT to quantify discrete- 
ness effect s in the usual case of a simple cubic (SC) lattice 
as preIC, iMarcod (|2008l ) develops fully its generalisation to 
the cases that preIC is a body centred cubic (BCC) or face 
centred cubic (FCC) lattice. We will exploit fully this latter 
generalisation in the next section. 

The principle of this approach is very simple: it con- 
sists simply in Taylor expanding the force on each particle 
due to any other in their relative (vector) displacement from 
the lattice configuratiorF^. Since the force is zero in the un- 
perturbed lattice the force F(R) on a particle originally at 
lattice site R can be written, at linear order in the displace- 
ments u(R, t), as 

F(R) = -^D(R-R')u(R',t), (3) 

R' 

where the sum is over all the lattice sites, and the matrix T> 
is 

D,.(R/0) = Gm(|^-3^) 

2?,-(0) = -^©^.(R) (4) 

where (5^,^ is the Kronecker delta, and the subscripts are the 
cartesian indicefl. With this approximation to the force, 
the equations of motion for the particles Eq. ([T]) may then 
be written as 

u(R, t) + 2Hu(R., t)^--^Yl ~ R')u(R', t) . (5) 

R' 

Deflning the discrete Fourier transform on the lattice and 
its inverse by 

Q(k, t)^J2 e"* '^u(R-, t) (6a) 

R 

u(R,t) = l^e*-^u(k,i), (6b) 

k 

where the sum in Eq. (|6b[l is over the first Brillouin zone 
(FBZ) of the lattice, i.e., the set of A'^ non-equivalent recip- 
rocal lattice vectors closest to the origin k = cEl, Eq. (O 

The treatment is ana logous to one used standardly in solid 
state physics (see e.g. iPin cs (1963)) to treat perturbations about 
a crystal, both for the case of short range two-body i nteractions 

I e.g. L ennard-Jones) and Coulomb interactions. See [Marcos et aP 
20061) for further discussion. 

^"^ A sum over the copies, due to the periodic boundary condi- 
tions, is left implicit in these expressions. 

18 For a SC lattice the vectors of the FBZ are thus k = n.(2n/L), 
where n is a vector of integers of which each component rii (i = 
1,2,3) takes all inte g er va lues in the range 

See iMarcosI l|2008h for the explicit expressions for the 
FBZ vectors of a FCC and BCC lattice. 



can be written in reciprocal space as 

u(k,^)-|-2i^(^)ii(k,^) = --i^©(k)u(k,t) (7) 

where 0(k), the Fourier transform (FT) of V{R), is a sym- 
metric 3x3 matrix for each k. 

The solution of the dynamical problem now reduces sim- 
ply to the diagonalisation of the ©(k), which is straighfor- 
ward (and inexpensive) numerically.. For each k this gives 
three orthonormal eigenvectors e„(k) and their eigenvalues 
a;^(k) (n = 1, 2, 3). The evolved displacements from any ini- 
tial perturbed lattice configuration, specified at a time to, 
may then be written as 

u(R, t) = ^ E [^(k, t)G(k, to) + Q(k, t)a(k, to)] e''''^ 

k 

(8) 

where the matrix elements of the "evolution operator" V 
and Q are 

3 

Pp4k,t) =E(7„(k,t)(e„(k))^(e„(k)), (9a) 

n=l 
3 

Q^4k,t) =EK(k,t)(e„(k))^(e„(k)),. (9b) 

n=l 

The functions C/„(k,t) and V„(k,t) are linearly independent 
solutions of the mode equations 

/ + 2ir/ = -^^/ (10) 

chosen such that L''„(k,to) ~ 1, f/n(k, to) = 0, Ki(k,to) — 
and V;i(k, to) = 1. 

The expression Eq. ^ for the evolution is, up to the 
validity of the linearized approximation to the force, exact 
for the discrete system. To use it to determine discreteness 
effects we must, as we have discussed at length in the pre- 
vious subsections, first identify unambiguously the correct 
continuum limit, and how it is obtained by extr apolation of 
the discreteness parameters. We have shown in lJovce et al.l 
(2005); Marcos et al. (2Q0^) that this may be done straigh- 
forwardly directly from Eq. ([8}: taking the limit ^ — > 0, 
this expre ssion for the e volution converges exactly to that 
obtained l|Buchertlll993 ) by linearizing the equations for a 
self-gravitating fluid in the Lagrangian formalism, which re- 
duces asymptotically to the Zeldovich approximation. The 
latter represents the appropriate analogous treatment of the 
Vlasov-Poisson limit. Note again that the limit ^ — > in 
Eq. ((SI is taken at fixed to (i.e. flxed initial red-shift Zi) 
and for a fixed input spectrum of displacements u(k) (and 
velocities u(k))Ej. 

The differences between the evolution given by Eq. ([SJ 
and this continuum evolution may then be computed exactly 

1^ The convergence to the VP limit is shown most easily keeping 
the initial conditions fixed, as described above, i.e., by introducing 
an ultraviolet cut-off in the initial spectrum so that the initial 
modes remains the same as £ — > 0. Alternatively, and in line with 
the general prescription given above, the limit may be recovered 
with a finite smoothing e which is kept fixed as £ — > 0. To do 
so one exploits the fact that the PLT formalism can be applied 
to any two body interactio n potential, and specifica lly a softened 
gravitational potential (see I Jovce fc MarcosI l|20073) for details). 
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for any given initial conditions, yielding the discreteness ef- 
fects, in the r e gime of validity of the PLT approximation. In 
iMarcos et alj (200(f) we have studied the domain of validity 
of PLT numerically for a range of different initial condi- 
tions, and conclude that, for some simple statistical quanti- 
ties, it describes the evolution up to the time when many 
particles approach their nearest neighbours, which corre- 
sponds approxi mately to shell cr o ssing i n the corresponding 
fluid limi1r°l. In ljovce fc Marco^ (j2007al ) we have presented 
a quantitative analysis of the discreteness effects in this cor- 
responding regime, for some basic quantities in typical cos- 
mological simulations. Some of the essential results are the 
following: 

• The modification at shell crossing (up to which the PLT 
treatment described fully the discreteness effects) of the evo- 
lution of any given mode of the displacement field grows 
monotonically with time. Indeed taking the limit 2^ — > cxd 
at fixed particle density (i.e fixed £) the evolution of the 
N-body system diverges from the continuum VP limit. 

• The modification due to discreteness at shell crossing, 
for a fixed Zi, depends approximately on the ratio k£, in- 
creasing as k£ does. This is physically very reasonable: the 
longer the wavelength of a mode compared to the interpar- 
ticle scale, the less affected is the evolution by discreteness. 
For the typical values of Zi in cosmological simulations, the 
effect is typically to reduce the power in modes, by up to 
about 50% at the Nyquist frequency, and by about 10% at 
half this value. 

Given this treatment of discreteness effects — exhaus- 
tive and analytical, but with a limited domain of validity in 
time — what can we conclude about the questions raised 
in the previous subsections? Concerning the formal extrap- 
olation to the VP limit, we have already noted that PLT 
indeed converges to the theoretical VP behaviour when the 
parameters are extrapolated as prescribed above. With re- 
spect to the question of how far we must extrapolate in £ in 
order to converge with some required precision to the VP 
evolution, the treatment also gives a clear answer, at shell 
crossing. The answer depends of course on the quantity con- 
sidered, and then also on Zi and on the cosmological model 
(which de termined the red-shift of shell crossing given £). In 
Ijovce &:~ Marcos (2007a) we have shown, for example, that, 
for Zi a factor of five larger than the redshift of shell cross- 
ing, errors of five percent in the PS are achieved only for 
k < fcjv/4. Thus if we want, at shell crossing, an accuracy of 
less than this on the power, we can use only results in this 
range. 

This second conclusion is strikingly different from what 
one might expect given the "common wisdom": the evolu- 
tion of the simulation makes the range of scale over which 
the continuum model is accurately represented (to some 
given precision) decrease, rather than increase. It in fact 
suggests that the view that £ should be a lower limit for 
spatial resolution may even be too optimistic. These find- 



We will see below that, for the quantities which we will study 
in the next section — the differences in the evolution of the same 
initial conditions sampled on different preIC — the PLT approx- 
imation actually hold for much longer times, apparently following 
well the evolution of any mode until it goes non-linear. 



ings, however, only apply at shell crossing, and the "com- 
mon wisdom" above may still apply later on. Indeed, as we 
described, the justification for this common wisdom is that 
when the transfer of power, characteristic of gravitational 
clustering in these systems, sets in, differences at smaller 
scales are wiped out. These results at shell-crossing show, 
however, that between Zi and shell crossing errors develop 
in the long- wavelength modes (below fe]v) which were not 
present in the initial conditions. As a result modes at later 
time which depend on this power will necessarily inherit this 
error. 

An important point which we emphasize is that the fun- 
damental reason why the discreteness errors determined us- 
ing PLT do not behave as expected by the common wis- 
dom is that they arise from physical effects of discreteness 
which are not usually envisaged. Indeed the physical effects 
described by PLT compared with the VP limit are different 
from the two effects envisaged usually which we listed above. 
Firstly, they are dynamical effects which modify the evolu- 
tion of any given mode in a way which is independent of the 
initial conditions. Secondly, they are clearly not two body 
coUisional effects 0- The effect they describe can be charac- 
terised physically as a dynamical sparse sampling effect: PLT 
compared with its VP limit tells us how the evolution of a 
fluctuation depends on the spatial density of the sampling 
particles. An important question is then evidently to under- 
stand how this physical effect — which there is no reason 
to believe should go away when we pass to the non-linear 
regime — quantitatively affects results in the latter regime. 
We will return to this point in our conclusions. 



3 A CALIBRATED NUMERICAL STUDY OF 
DISCRETENESS EFFECTS 

We return now to the practical question of how small £ needs 
to be for a measured quantity to have converged to a desired 
precision. Since the force smoothing e places a lower bound 
on the spatial resolution, a simplified, more specific, form of 
the question is: how small does £ have to be in order that, 
at any given red-shift, the effects of discreteness are negligi- 
ble down to scales of order e? The answer provided by the 
"common wisdom" above is that £ is sufficiently small, in 
typical simulations, i f £/e i s less than about one hundred 
(see e.g. iKnebe et al.l l|2000t )). According to the "dissenting 
views" £ must be at least as small as e. 

One way to determine, in principle, which view is cor- 
rect is evidently to compare results from simulations with 
large £i/e in the range £i > r > e with those obtained in 
much higher resolution simulations, with £2 ^ e <^ £ i. Thi s 
is indeed the strategy advocated in ISplinter et al.l l|l998l ). 
which reports a study of this type down to a resolution 
£2 = e. It concludes, as noted above, that there are signifi- 
cant differences in results, i.e., no eviden ce for co nvergence, 
in the range £1 > r > £2. Other authors (|Knebe e t al. 200_3) 
argue however that these differences are ascribable to "er- 
roneous evolution in high resolution runs". The difficulty 

To m ake this very explicit we have shown in ljovce &: MarcosI 
ll2007all that the inclusion of a simple Plummer smoothing in the 
force actually increases the difference between PLT and the VP 
limit for unsmoothed gravity. 
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in reaching a convincing conclusion is that the questions of 
discreteness effects are intertwined with numerical and finite 
size effects. While such differences should be resolvable by 
further numerical tests, this would require considerable in- 
vestment of resources which, apparently because of the wide 
acceptance of the "common wisdom" , has not been mad43- 

Instead of undertaking such a numerical study — which, 
given the modest numerical resources at our disposition, 
would not in any case likely to b e any more conclusive than 
that reported by ISplinter et al.l l|l998t ) — we focus in the 
rest of this paper on another kind of test. We will see that 
this will allow us to reach conclusions, with modest sized 
(but very well numerically converged) simulations, about the 
central issue: the validity /precision of results in the range of 
scales around or below I, in simulations with e. The aim 
is to provide a method which gives a non-trivial lower bound 
on discreteness error in such simulations. To do so we simply 
compare the results of simulations from identical theoretical 
initial conditions, changing only the choice of the discrete- 
ness parameter preIC, i.e., the pre-initial configuration. We 
can then study how this error depends on time and scale. 
Although the measured effects are quite small — at most of 
order of five percent in the PS for the times and scales rele- 
vant to cosmological simulations — we can establish clearly, 
using the analytical PLT formalism combined with numeri- 
cal tests of their dependence on £ and e, that they are indeed 
discreteness effects. We can then address in a controlled way 
the question of how far £ needs to be extrapolated so that 
one can be confident that the true systematic errors due to 
discreteness have converged to significantly less than this 
lower bound (e.g. to less than one per cent). 

Rather than considering a specific cosmological model, 
we consider a simple power law PS with exponent n=-2, 
evolved in an EdS universe. This choice is both suitable for 
our study as it is simple — introducing no characteristic 
scale in the input model — and yet close to the currently 
favoured CDM-like cosmological model, which has an ini- 
tial PS with effective exponent ranging between n ~ — 1 
and n « —3 over the relevant range of scales. In particu- 
lar we note that this PS is, like these cosmological models, 
long-wavelength dominated so that the very efficient trans- 
fer of power from long to short wavelengths which, as we 
have discussed above, is believed to play a role in wiping 
out discreteness effects, should be well represented. We will 
comment further in our conclusions on the generalization to 
other initial conditions, and specifically to those of currently 
favoured cosmological models. 

All our simulations have been performed using the 
publically available parallel tree-mesh code GADGET2 
jSpringel et al.ll200ll ). We use this single (widely used and 
highly tested) code for our study for the reasons we discussed 
above: the discreteness effects we are trying to understand 
and control for are distinct from differences arising between 
different codes, and indeed distinct from any dependence of 
results on the numerical parameters of a given code. The 
"calibration" of our results with our analytic tools here pro- 
vide in fact a robust check that the GADGET2 code's in- 
tegration of the iV-body equations of motion is sufficiently 



Figure 1. From left to right, unit cell of the SC, BCC and FCC 
lattices. 



PI configuration 




N 


SC 




= 262144 


BCC 


2 X 51^ 


= 265302 


FCC 


4 X 40^ 


= 256000 


glass 


64^ 


= 262144 



Table 1. Number of particles in the four PreIC of the reference 
set of simulations SI 



precise that this is indeed the case. Comparison with other 
codes would be, in the relevant regime, a check on the ac- 
curacy of these codes, rather than a check on our results. 
In the regime where our analytic results do not apply, we 
can have, of course, less confidence in the identification of 
our measured effects as physical discreteness effects, and a 
comparison with other codes could be instructive. We will 
address this issue below, where we give details of the de- 
tailed checks of numerical convergence of our results which 
we have performed using GADGET2. 

3.1 Initial conditions 

We use the standard method, based on the Zeldovich ap- 
proximation, to set up initial conditions by applying appro- 
priate displacements to four different preIC: a simple cubic 
(SC) lattice, a body centred cubic (BCC) lattice, a face cen- 
tred cubic (FCC) lattice, and a glass configuration, shown 
in Fig. [U 

Our reference set of simulations, which we denote SI, 
have the number of particles shown in Table [T] The num- 
bers for the BCC and FCC configurations have been cho- 
sen to be as close as possible to those of the SC and glass 
configuration^^- The glass is generated, starting from Pois- 
son distributed points with zero velocity, using an option in 
the GADGET2 package which evolves the particles under 
Coulomb forces (without expansion) and with a damping 
implemented by setting the velocities to zero at each time 
step. In what follows our results will always be given in units 
of length m which the box size is equal to unity. In these units 
the value of £ in the four different preIC varies by less than 
one percent. 

The three lattices have PS which can be written 



P(k) = (2^) 



K) 



(11) 



See, however, the recent paper bv lRomeo et al.l ll2008a^ . which 
we will comment on in our conclusions. 



As can be seen in Fig.[T] there are two particles per elementary 
cell of a BCC lattice, four per cell in a FCC lattice. Thus in a 
cubic box we have 2M^ in a BCC, 4M^ in an FCC, lattice, M is 
an integer. 



10 M. Joyce, B. Marcos and T. Baertschiger 




Figure 2. The power spectmm(PS) of the pre-initial glass. The 
dashed line indicates the behaviour oc fc'*. 



where S(k) denotes the Dirac delta function and the sum 
runs over non-zero vectors K which are an infinite subset of 
the full reciprocal lattice appropriately defined for a given 
lattic43- The delta function structure of these PS is a result 
of the translational symmetries of the lattices. The PS of 
the glass is, in contrast, a continuous function. Indeed, up 
to finite size effects, it is a function only of fc = |k| because of 
its statistical isotropy. It is shown in Fig. O along with a line 
indicating its approximate small k behaviour, P{k) oc A;"*. 

Given an input theoretical PS Pth{k), we generate a 
realisation of the displacement field u(x) to be applied to 
the particles at spatial positions x of the four preIC in a 
cubic box taking 

u(x) = [ak sin(k ■ x)k — 6k cos(k ■ x)k] , (12) 



with 



flk 



R 



where _Ri and R2 are two independent Gaussian random 
numbers with dispersion equal to unity. In writing the dis- 
placement field as a Fourier sum we use the fact that the 
preIC are set up on a periodic cube, and the sum over the 
vectors k extends then over the appropriate reciprocal lat- 
tice. Further, if the input PS itself does not have a cut- 
off at a wavenumber significantly smaller than the Nyquist 
wavenumber of the sampling preIC distribution, such a cut- 
off must be imposed to avoid aliasing effects. Here, where we 
consider a simple power-law PS without a cut-off, we will 
take the sum in k to extend over the first Brillouin zone of 
the SC lattice, i.e., the reciprocal vectors k = n(27r/L) with 
each integer component rii £ [—N/2,N/2[. As we will dis- 
cuss further below, this is the choice which minimises alias- 
ing effects for the SC lattice, but not for the other preIC 
configurations. We will measure the associated very small 
aliasing effects in the initial conditions and keep track of 



2'* For the SC lattice we have simpl y K = n(27r/£ 3c) where n is 
any non-zero vector of integers; see iMarcosI ll200d ) for the more 
general definition for any Bravais lattice. Note that when we use 
the term "reciprocal lattice" here, we are refering to that defined 
for the periodic box of side L, i.e., for the SC lattice K = n(27r/L) 
where n is any vector of integers. 



their role in generating differences in the evolved distribu- 
tions. 

The only other parameter which needs to be fixed is the 
normalization of the input PS Pth(k) (which is equivalent to 
the choice of the initial red-shift Zi). In the set SI we have 
taken, for all preIC, 



k%P{kN) = 0.6, 



(14) 



where, in our units, fcjv = 647r. We have made this choice for 
our reference simulations because it is close to that chosen 
for such initial conditions, on a SC lattice, by the widely used 
GRAPHICS package (Bertschinger 19950 We wiU discuss 
below the effect of modifying this choice. 

Before turning to the evolution from this set of initial 
conditions, let us consider more precisely their correlation 
properties, and in particular the effects of aliasing we have 
mentioned above. To do so we make use of the detailed anal- 
ysis of initial co n ditions of N-body simulations reported in 
iJovce fc Marcos! (l2007bf ). The PS of the perturbed preIC 
distribution can be written conveniently in the form 



P(k) =Pc(fc)+Pd(k). 



(15) 



where Pdk) is the "continuous" part, independent of the 
preIC distribution, and Pd(k) is a "discrete" term which 
depends on the latter. The full analytic expression for both 
these quantities can be expanded order by order in the am- 
plitude of the input theoretical PS Pth(fc). At leading order 
one obtains 



J'c(fc) 

Pd(fe) 



Pth(fc), 

Ppl(k) + P,l(k) 



(16) 



where Ppi(k) is the PS of the unperturbed preIC configu- 
ration (i.e. lattice or glass) and 



P.,(k) = ^ / '^'?(q-k)'^^[Ppi(q + k)-Ppi(k)] 

(17) 

is a contribution to the PS which, if non-zero at small fc, de- 
scribes an aliasing of the input PS. If the preIC is a perfect 
lattice, Ppi(k) is given by Eq. ((TTJ and therefore, for k ^ K, 
we have 



[(K - k) ■ k] 

P.i(k)=fc^^ L J Ph(K-k), (18) 



where K are the appropriate subset of reciprocal lattice 
vectors for each lattice. For the SC lattice the vectors K 
are given by K = = 2fcivn, where n is a ny non-zero 
vector of integers. It can then be verified easiljifj from the 
expression in Eq. (|18p that Pai(k) is zero inside the first 
Brillouin zone of this lattice, i.e., for all k with each com- 
ponent ki g] — kisi,kjsi], if we impose a cut-off by making 
Pth(k) zero outside the same region. This choice, which is 
the one we have used (and which is that standardly used in 

25 This package determines the starting red-shift Zi by normal- 
izing so that the maximal value of the density fluctuation field 
at any point of the lattice is unity. This gives a mass variance 
at this scale considerably less than unity, sufiiciently small that 
non-linear corrections to the Zeldovich approximation should be 

small. 

° See Ijovce &: Marco i ll2007al) for a more detailed discussion. 
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this context lBertschinge3 (|l995l ): ICouchmanI (|l99ll ')') IS op- 
timal, in the sense that it maximizes the size of the region 
about k = in reciprocal place where the representation of 
the input power is exact, at linear order in the amplitude of 
the input PS. 

For other lattices an analagous, but different, optimal 
choice can be made, taking the input PS non-zero only inside 
the given lattice's first Brillouin zone. Here we have not done 
so as such a procedure would require sampling the input PS 
at different wavevectors, which is incompatible with the re- 
quirement that we use an identical realization of the theoret- 
ical initial conditions. As a result we will have for the non-sc 
lattices a small contribution coming from the aliasing term 
as given in Eq. (|f 8p . For the glass, on the other hand, a more 
significant contribution from this term [as given in Eq. (I17p ] 
is expected, as it is alwa ys non-zero a nd proportional to at 
sufficiently small k (see .Joyce fc Ma rcos (2007b) for further 
detail). In what follows we will study carefully in simula- 
tions the evolution of these residual differences at small k 
power in the initial conditions, showing that they can in fact 
be neglected in understanding the differences in the evolved 
power which emerge at these scales. 

3.2 Numerical Evolution of SI 

We evolvj^ these four initial conditions in an EdS cosmol- 
ogy, from a scale factor a = I to a = 2^. At this final time, 
as we will see below, the scale of non-linearity has reached 
the box size and finite size effects dominate. GADGET im- 
plements a smoothing which modifies the force from exactly 
Newtonian only below a scale e. We take here e = ^sc/15, 
where ^ac is the interparticle spacing of the SC lattice. This 
is, according to the "common wisdom" , a conservative choice 
for the final resolution scalCJ. 



3.3 Snapshot inspection 

In Figs [3] and |31 we show snapshots, for each of the four 
initial conditions, of a slice of depth 0.31/ of the simulation 
box. The four snapshots correspond to a = I, a = 2^, a = 2^ 
and a = 2^. In the initial conditions, at a = I, the dis- 
tributions look very different, reflecting the different small 
scale properties, and long-range order, of the preIC con- 
figurations. Blurring slightly one's vision, however, one can 
make out clearly in the lattice configurations the very similar 
superimposed fluctuations at larger scales. The glass looks 
very different because it does not have the deterministic long 
range order of the lattice, which makes the projection ap- 
pear considerably densei|fj. In the second slice, at a = 2"^, 

The details on the numerical parameters we have used for the 
results reported are given in Appendix 1X1 

For comparison we note that, if our comoving particle den- 
sity is assumed equ al to that in the Millenium simulation 
llSprineel et al.l 1I2OO5I) '). the comoving size of our box is then ap - 
proximately 15 h ~^ Mpc. The ratio £/e in ISpringel et al] tOok) 
is approximately fifty. 

In passi ng we underline that, contrary to what is sometimes 
stated (e.g. IWang fc Whit j \200l\) ). the glass is a long-range or- 
dered distribution. In fact it has the property that P{k = 0) = 0, 
which imposes the global constraint that the integral of the two 
point correlation function is zero. Discussion of the very partic- 



the flrst non-linear structures have formed, and already now 
the visual impression is of a very strong resemblance in the 
clustering. Distinct differences are however still evident. In 
particular alignments inherited from the lattice conflgura- 
tions are clearly visible, most evidently in the SC lattice. 
In the next slice at a = 2® (which, as we will see below, is 
about the time at which the largest modes included in the 
box go non-linear) the first visual impression is of an even 
greater resemblance of the configurations, but again closer 
inspection reveals differences at smaller scales. Likewise in 
the last slice, when almost all the mass is in just a few halos, 
the broad features at large scales are impressively similar, 
while the spatial organisation of smaller structures reveals 
evident differences. 

Some of the differences observed visually in the earlier 
time snapshots are manifestly related to the subtle differ- 
ences in the initial conditions, and are therefore clearly dis- 
creteness effects. The differences in the more evolved snap- 
shots are, however, not necessarily indicative of anything 
other than the intrinsically chaotic dynamics of the non- 
linear regime of the evolutiorP°l. What we are interested in, 
and will now examine, are differences in the statistical prop- 
erties of these distributions, which are what we use them to 
infer in cosmology. 

3.4 Power spectrum and correlation function 

Let us consider more quantitatively the differences in the 
two point properties of these distributions. In Figs. [5] and [6] 
we plot the reduced two point correlation function £(r) and 
the PS P{k), for a series of four different time sliceijfj. Also 
shown, in an inset panel in each case, are the normalised 
residuals of each quantity with respect to the average, i.e., 
for a quantity C^(i) in the i-th bin (of k or r) in the simu- 
lation of initial conditions I {I—SC, BCC, FCC, glass) the 
residual is 

^c'{i)-j:,c'{i) 



5C'{i) 



(19) 



Also shown in each case is the "linear theory" (LT) predic- 
tion for the evolution of the theoretical PS, i.e., the initial 
theoretical PS multiplied by . 

These results reflect broadly the impression gained by 
visual inspection above. In particular inspection of the cor- 
relation function shows four distributions which are appar- 
ently very different at the initial time evolve to closely resem- 
ble one another already at a = 2^. Note indeed that, at this 



ular stochastic long-range order of such "superhomogeneous" (or 
" hyperuniform" ) distributions , with P(k = 0) = 0, may be fo und 
in lGabrielli et al.l 1I2OO2I . I2OO3I ): iTorguato fc StillingeJ ll2003t) . To 
generate them starting from a Poisson distribution, as here, one 
requires long-range correlation in the displacements of the parti- 
cles, provided here by the dynamics under Coulomb force (which 
rearranges the points so that the fluctuations in any volume are 
proportional to the surface i.e. sub-Poissonian). 
^'^ For a discussion of chaos in_/V-body self-gravitating systems 
see, e.g.. lSideris fc Kandrupl l|2002h . 

Details of how these quantities have been estimated are given 
in Appendix |B] Note that we write the two point correlation 
function and the PS as functions only of the modulus of their ar- 
guments as an average over spherical shells in real and reciprocal 
space respectively is performed. 
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Figure 5. PS (left column) and correlation function in real spa<;e (right column) for a = 1 (top row) and a = 2^ (bottom row) 



time and the subsequent ones, the correlation functions are 
so similar as to be indistinguishable in the main plot down to 
about r = 10~^, which coincides with our chosen e. In recip- 
rocal space the resemblance of the initial conditions — the 
fact that they represent exactly the same realisation of the 
input PS — can be seen. Indeed the initial power below fcjv 
agrees in all cases to a precision of less than a small fraction 
of a percent. Already in the next time slice shown, at o = 2^, 
the PS in the main plot of the four distributions are super- 
imposed almost perfectly over the entire range, except for 
a still visible difference for the sc configuration, correspond- 
ing to the difference we identified by visual inspection. In the 
last time slices the curves for the PS are again, as in the case 
of C(r), effectively indistinguishable over, in the final slice, 
more than four orders of magnitude in power. Note that 
in our length units the asymptotic Poissoniam behaviour (of 
any translationally invariant stochastic point process) corre- 



sponds to P(k ^ oo) = 1/no = 1/64^ ^ 4x 10"^. When the 
PS asymptotes to this value it indicates that the PS mea- 
sured is dominated by the intrinsic noise of the discrete pro- 
cess. We see that this maximal resolved wavenumber propa- 
gates to larger k in time, reaching a final value of order n/e 
(i.e. k/kN = £/e). 

This first view of these results thus supports very 
strongly the "common wisdom" : there is very efficient trans- 
fer of power from large to small scales which wipes out 
memory of the differences between the initial conditions at 
small scales. Rapidly results converge down to a scale char- 
acterised by e. There is no significant dependence on preIC, 
and therefore the associated discreteness effects are wiped 
out. 

A more careful analysis of the evolution of the resid- 
uals in each plot shows, however, some behaviours which 
are unexpected according to this common wisdom. Firstly, 
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Figure 6. PS (left column) and correlation function in real space (right column) for a = 2^ (top row) and a = 2^ (bottom row). 



for fc < fcjv we see differences in the PS which appear and 
grow monotonically with a, i.e., the gravitational evolution 
appears to produce some small differences at large scales 
which were not present initially. While suc h an effect is pre- 
dicted by PLT, as documented in detail in I Joyce fc Marcoa 
l|2007al ) and described briefly above, this is valid only at suf- 
ficiently early timeiP^. Thus the expectation that non-linear 
transfer of power from larger scales may wipe out the effects 
of PLT at these times appears not to be correct. For k > kN, 
on the other hand, such an effect is indeed apparent, but ap- 
pears only to be operative at the very earliest times when 
the non-linear structures first develop at smaller scales. In- 
deed, between a = 2^ and the final plot, at a = 2^, there 



Note that since the residuals are normalised, they are constant 
under fluid LT. 



appears to be little evidence for any further washing out of 
the residuals in the power. On the contrary they appear to 
grow, giving a dispersion of order several percent at the final 
time. 

These results are clearly illustrated in Fig. [T] which 
shows the square root of the variance of PS calculated at 
each fc, and each of the four times as indicated, over the 
four realizations, i.e.. 



api{k,a) 



1/2 



(20) 



where I labels the m = 4 different preIC, and we have 
defined 



(21) 
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Figure 7. Normalized variance of the PS <Tpj{k), defined in 
Eq. |[20t . 



Given that the amplitude of the effects are so small — at 
small k and early times in particular — we evidently need 
to be careful in interpreting these differences as resulting 
from the physical discreteness effect we set out to measure. 
We could envisage that such a time and space dependent 
dispersion could be the produced, in particular, by numer- 
ical effects in the evolution or by statistical effects in the 
estimators. For example, it is conceivable that there is an 
interplay between the numerical errors relating to the cal- 
culation of the force and each particular initial condition, or 
that the variance measured is simply a statistical variance 
which would decrease if we took more particles (i.e. a larger 
box at the same particle density). In the rest of this section 
we examine this question carefully, establishing — we be- 
lieve very convincingly — that, at least up to the slice at 
a — 2^ , this measured dispersion is a discreteness effect. 



3.5 Numerical convergence 

Let us first consider the stability of the results with respect 
to variation of the numerical parameters, i.e., those control- 
ling the accuracy of the numerical integration of the A''-body 
equations at given values of the discreteness parameters . 
In the GADGET2 A'^-body code, there are two sets of such 
parameters: a first set controlling the time-stepping and a 
second one the resolution in the calculation of the force. In 
Appendix |X] we give the full details of two sets of param- 
eter choices for which we now compare results: a "low res- 
olution" (LR) simulation, corresponding to the values used 
in obtaining the results given above and subsequently in 
the paper, and a "high resolution" (HR) simulation. As we 
discuss in further details in Appendix [Xj the LR are typi- 
cal choices for large cosmological simulations in the litera- 
ture (e.g. tho s e of t he VIRGO consortium, as described in 
IJenkins et all lIlQQSl )). while our HR values are even more 
stringent choices than typically use d in similar converge nce 
tests reported in the literature (e.g. ICrocce et ah! l|2006f )). 

We illustrate the degree of convergence between the LR 
and HR simulations in Figs. [8] and |9l The former shows 
the excellent stability of normalized differences like those 
we have considered above — and will focus on in the rest 
of the paper — in the PS for the two preIC indicated. The 
latter figure shows, on the other hand, the differences of the 
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Figure 8. Normalized diff-erences of the PS for the BCC and SC 
preIC configuration, for high resolution simulations (thick lines) 
and low resolution simulations (thin lines) at a = 2^, a = 2^ and 
a = 2'^. 



results of the LR and HR simulations for the full PS in each 
of two preIC taken separately. We see that these differences 
are, at the two later times, comparable in magnitude to the 
differences we measure (in the previous figure) , over a part of 
the range of k. Thus the full PS measured in each of the two 
preIC simulations changes as a function of the numerical pa- 
rameters in this range by as much as the differences between 
them which we are studying here (and which we have just 
seen to be well converged numerically). This means simply 
that the numerical errors associated with these changes in 
parameters are correlated strongly with the full PS, which 
is very close to the same in the two cases, and so cancel 
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Figure 9. Normalized differences of the PS of HR and LR runs, 
at a = 2^, a = 2^ and a = 2^. In each panel the thick lines 
corresponds to the BCC prelC configuration and the thin one to 
the SC one. 



out when we take the difFerence. This suggests that, in more 
general, it may be easier to place this kind of lower bound 
on discreteness effects than to attain a comparable level of 
numerical convergence on other quantities (such as the full 
PS). 



3.6 Comparison with PLT 

The PLT formalism for the evolution of the displa c ement s 
off the lattice , deve loped explicitly in Ijovce et al.l (|2005h : 
iMarcos et al.l (20061) for the SC lattice, has been general- 
ized in iMarcosI l|2008l 'l to both BCC and FCC lattices. We 



exploit these analytical results here, for the case of the SC 
and BCC lattice, as a control on the accuracy of our nu- 
merical simulations at sufficiently early times when PLT is 
a valid approximation. Conversely this comparison can be 
seen — given the results just shown above on the numeri- 
cal convergence of our results — as a check on the range of 
applicability of PLT. We will see that this range turns out 
to be considerabl y greater than that which was established 
in the studies in [Marcos et al. I (|2006l ). making PLT a very 
useful tool for calibrating numerical results. 

To compare our numerical results with PLT we simply 
generate, for each set of BCC and SC initial conditions, the 
configurations given by PLT evolution of Eq. (IS}, where the 
eigenvalues and eigenvectors are those for the corresponding 
la ttice. The deta ils of these latter calculations may be found 
in lMarcod (|2008l ). 

In Fig. [To] are shown, for the GADGET2 simulations 
and the PLT evolved configurations, the normalised differ- 
ences between the PS for the BCC and SC, i.e.. 



Pbcc(fc, g) - Psc(fc, g) 
i [Pbcc(fe,g) +Psc(fc,g) 



(22) 



where the subscript indicates the prelC. We also show for 
comparison in Fig. [11] the same quantities, except that the 
PLT evolved configurations are replaced by those evolved 
with its fluid limit (which we will denote by FLT, for "fiuid 
linear theory"). For initial conditions set up, as done here, 
with the Zeldovich approximation, this is simply the extrap- 
olated evolution in this same approximation. 

The agreement with PLT at g = 2'^ is extremely good 
for all the measured fc, while at g = 2^ it is restricted only to 
the very longest wavelength modes in the box. FLT, on the 
other hand, traces the observed differences well until g = 2'', 
but only the k larger than fcjv. 

These different ranges of agreement for PLT and FLT 
are simple to understand, using the results quoted above 
in Eqs. (|15ll7p . These formulae relate, at sufficiently early 
times and small k, the theoretical PS of density fiuctuations 
fth(fc) to the full PS of density fiuctuations in the generated 
point distributions. FLT gives a linear amplification of the 
displacement fields, independent of k, and therefore a linear 
amplification of the terms -Pc(k) and Pai(k). Outside the 
range of fc where Ppi(k) contributes, i.e., inside the FBZ, 
FLT thus simply describes a linear amplification of the full 
initial PS, which leaves the normalized quantity in Eq. (|22|l 
strictly invariant. Outside the FBZ, on the other hand, the 
term Ppi(k) becomes important. When this is the case the 
full evolution is well approximated b y the FLT evolution 
because (see jjovce fc Marcos! (|2007al )) the evolving term, 
Pai(k), is in fact dominated by initial power at small k for 
which the evolution is very well approximated by FLT. 

The regime in which PLT traces the differences very 
well, but FLT does not, corresponds to the k inside the FBZ 
which are, in PLT, amplified linearly in slightly different 
ways on each lattice. In this case the physical discreteness 
effect arises thus from the modification with respect to FLT 
of the dynamical evolution of the same initial power. In the 
regime in which FLT gives a good approximation, on the 
other hand, the corresponding discreteness effects arise from 
the power associated with the slightly different initial sam- 
plings on the different lattices of modes which evolve ap- 
proximately in the same way. 
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Figure 10. Normalized variance of the PS of the BCC and SC, computed from our GADGET simulations of full gravity (continuous 
lines) and simulations evolved using PLT (dashed lines) for, from top to bottom and left to right, a = 2", a = 2^, a = 2^ and a = 2^. All 
the curves are normalized to the PS for the full gravity (FG) case. 
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Figure 12. Normalized difference between the PS computed with 
FG and PLT for the BCC (thick lines) and SC (thin lines) con- 
figurations. 



Let us consider further the range of validity of PLT 
in these plots. The perturbati v e exp ansio n underl y ing P LT 
as developed in iMarcos et"aLl (|2006l ) and iMarcosI (|2008l ) is 
strictly valid, as we have discussed above, only when the 
relative separation o f all particles i s smal l compared to their 
initial separation. In lMarcos et all (|2006D it has been shown 
that it gives a very good approximation to the evolution 
of the PS (and significantly better than FLT) at least un- 



til the time when a significant fraction of the particles have 
come close to another particle for the first time (which cor- 
responds approximately to shell crossing in the fiuid limit). 
However, its possible validity beyond this time has not been 
established. What the results in these plots show is that its 
validity indeed extends considerably longer, as there has al- 
ready been very significant shell crossing already at a = 2'^, 
and clearly at a — 2^ the evolution is well beyond this point. 
In Fig. [T2] we show the modulus of the ratio 

Pr(fc,a)-Pf^^(fc,a) 

Pf°(fe,a) ^''^ 

for both I =SC and I =BCC, i.e., the fractional deviation 
of the power at each k in the PLT evolved initial conditions 
simulation from that in the full gravity (FG) simulation of 
the same initial conditions. Comparing with the results of 
Fig. IIUI we see that the range in which PLT correctly de- 
scribes the differences between the SC and BCC simulations 
extends in fact to when the plotted quantity is of order two, 
i.e., into a regime in which PLT no longer follows well the 
full power in each PS accurately. This is evidently possible 
only if the deviation from the PLT evolved initial conditions 
is essentially the same for both preIC, i.e., this additional 
non-linear power itself has smaller discreteness corrections 
than those given by PLT. We note that this is very coher- 
ent with our observations above concerning the numerical 
integration: we observed in that case much better numer- 
ical convergence of the measured differences than in each 
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Figure 11. Normalized variance of the PS of the BCC and SC PI, computed from FG simulations (continuous lines) and from simulations 
using FLT (dashed lines) for, from top to bottom and left to right, a = 2", a = 2^, a = 2^ and a = 2^. All the curves are normalized 
with the PS obtained with FG. 




Figure 13. Real deviation normalized by the estimated deviation 
computed using Eq. ((24}. 



of the PS individually. Thus these numerical residuals are 
strongly correlated with the non-linear power which is the 
same in both simulations, and so cancel out when we take 
the difference. We will discuss briefly in our conclusions these 
observations about the regime of validity of PLT. 

It is instructive also to examine, in the range in which 
PLT traces accurately the differences in the evolution, what 
the relation is between this measured difference, and the true 
discreteness error, which can also be calculated in PLT. In- 



deed, in the FBZ, it is simply the difference between the PLT 
evolved power and FLT evolved power. A simple qualitative 
measure is thus: 



^(fe)-n. 



Dev(fc, a) 



(24) 



w-pi, 



1/2 



where the a-dependence on the right hand side is left im- 
plicit. Limiting ourselves to the modes for which PLT fur- 
nishes a good approximation to the full evolution of the indi- 
vidual PS, i.e., to the regime in which the quantity plotted in 
Fig.[T2]is less than or order one, we see that Dev(fc, a) shows 
a clear tendency to increase with a, particularly for smaller 
k. Indeed at a = 2^, for the very smallest k for which PLT is 
still approximately valid, our lower bound on the discrete- 
ness error is one order of magnitude smaller than the real 
discreteness error. The reason is simply that the difference in 
the exponents characterising, in PLT, the growth of the dis- 
placement fields in these two different lattices at these values 
of k is considerably smaller than the difference between these 
exponents and the FLT behaviour (giving growth in propor- 
tion to a^). Even if, taking this factor into account, we arrive 
at a real discreteness error of order of one percent, this result 
shows that, for all our results here, we should bear in mind 
that the measured differences provide only lower bounds on 
the discreteness error which may be very different from the 
full discreteness error. 
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3.7 Variation of particle number TV 

A further direct numerical check on our interpretation of the 
differences we have identified as discreteness errors may be 
given by looking at their dependence on particle number A'' 
(i.e. effectively, given that we work at fixed box size, on the 
particle density parametrized by £). We thus consider vary- 
ing A'^ while keeping all other discreteness parameters fixed 
(and, again, checking also the stability of results considered 
to the variation of numerical parameters) . Shown in Fig. [TJ] 
are the results for the normalized variance api{k) [as defined 
in Eq. (|20|l ] on a set of four simulations of identical initial 
conditions (i.e. identical modes of the displacement fields 
on the four preIC) for iV = 32^ and N = 64^ particleQ 
Note that, to have identical displacement fields in the two 
cases, we have cut the initial PS at the Nyquist frequency 
of the iV = 32^ distributioiH- We see clearly the explicit N 
(or £) dependence of the results in all but the very strongly 
non-linear regime. In the PLT regime the difference in power 
depends parametrically on A'^ as N~^^^ [see Eq. (I26|l below]. 
Interestingly one can observe between a = 2^ and a = 2^ an 
apparent "spreading" of this explicit dependence to larger 
fc, a behaviour which is naturally interpreted as the transfer 
of the discreteness efi'ects accumulated at a = 2^ in the lin- 
ear regime to larger k modes as the corresponding scales go 
non-linear. On the other hand, we see no clear evidence for 
a dependence on TV at the strongly non-linear scales — and 
most notably over the entire range at the final time, a = 2^ 
— and so we will not assume here that the measured dif- 
ferences are discreteness effects. It is to be noted, however, 
that this is a very conservative assumption: the differences 
even at a = 2^^ may quite consistently be, and indeed are 
naturally, ascribed to those at the previous time, a — 2'^ , 
without the latter having to show the same explicit depen- 
dence on I. Indeed in the strongly non-linear regime |fj we 
do not expect a simple dependence of the final power on 
the amplitude of the power in the preceeding weakly non- 
linear phase, and therefore the dependence on I inherited 
from this phase could quite possibly be much weaker than 
that observed in the preceeding phase. 

3.8 Variation of e 

Another check on our results is given by considering the 
effect of varying e, keeping all the other discreteness param- 
eters fixed. In Fig. [15] we show again the normalized vari- 
ance api{k), now again for four different simulations with 
TV — 64'', for three different values of e: the same one as used 
in the results reported until now (e = £/15), and now also 
for simulations (from exactly the same initial conditions) 
with e = £ and e = 2£. We show only the range of k below 
the Nyquist frequency as this is the regime of physical inter- 
est, i.e., in which results are expected to converge to those 
for (unsmoothed) gravity, fixed approximately by the mode 

As described in Sect. 13.11 above, the number of particles in 
the non-SC configurations are chosen as close as possible to these 
numbers (Table [TJ. 

•^^ More precisely, as described in Sect. l3.l1 we sample exactly the 
modes in the FBZ of the SC lattice. 

We recall that, at a = 2*^, the whole box has gone non-linear 
with most of the matter in only a few halos (see Fig.|4j. 



inverse to the largest value of e. (We do not show results 
for smaller e as they are negligibly different in this range 
from those at e = £/15). The behaviour observed at a = 2^ 
is completely consistent with what is expected given that 
we have seen that PLT provides an excellent description of 
these differences at this time: the exponents for growth of 
the modes of the displacement field calculated in PLT (which 
may be calculated for any two-body potential) only begin to 
change significantly when e ^ £, simply because PLT is an 
expansion about the particles placed at their lattice sites. 
As e increases the deviation from the flu id evolution be- 
comes i n fact more and more significant (see lJovce fc Marcos! 
(|2007al )). but this deviation does not manifest itself as a dif- 
ference between evolution on the different lattices as the 
smallest scales on which they differ are then smoothed over. 
Thus the differences we measure decrease (in the FBZ, where 
they are due to the difference in the exponents relative to 
their FLT values). At a = 2^ we see essentially the same 
behaviour for the modes for which PLT was valid, while for 
the larger modes there is also some more marked decrease 
already for e = £. At a = 2^^ we see a larger spread, with an 
apparent tendency for the largest e to lead to the smallest 
differences, which would certainly be consistent with the hy- 
pothesis that these errors could also be interpreted as due to 
discreteness. It is important to note that, in all these figures, 
the reduction of the differences measured as e is increased 
does not imply a convergence of the simulations towards the 
physical (VP) limit, but at most towards a smoothed ver- 
sion of it, which may be further from the physical limit than 
the results obtained with the smallest e. Indeed in the PLT 
regime we have shown in IJovce fc MarcosI (2007a) that in- 
creasing e at fixed £ does indeed increase the deviation of 
the growth exponents of modes from their fluid value. 

3.9 Variation of initial red-shift 

The initial red-shift Zi is the remaining parameter in the list 
of discreteness parameters T^a we gave in our discussion in 
Sect. [5] As the dependence of varying it while keeping the 
other parameters fixed can be understood analytically using 
PLT, in the regime in which we know it to be valid (of small 
relative displacements), we do not report here numerical re- 
sult^3- Quite simply we note that, in the E dS cosmology, 
the evolution of the PS in PLT can be written l|Marcos et al.l 
(,2006 )') to a very good approximation as 

P(k, a) = a'+^"<'')*='' V(k, a = 1) (25) 

where 5i(k) is a function of the orientation of the vector 
which depends on the prelC. It follows that the normalized 
difference in the power, averaged in a bin of wavevectors 
centred at wavenumber k, scales approximately as 

~ [4c(fc)-5bcc(fc)] fc'^'^loga (26) 

where 5i{k) are appropriate effective values of the param- 
eter 5i(k) over the bins of wavevectors. The differences we 
have measured thus increase without limit as Zi does, with 
a logarithmic depedence on the latter. 

See also Ijovce fc Marco^ l l2007al ) for quantitative results. 
So me numerical results fo r the effect of varying Zi only are given 
in [McDonald et ah! ll2006h , but only for very specific quantities 
(ratios of PS for different dark energy models). 
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Figure 14. Normalized variance of the PS for TV = 32^ and = 64^ particles at, from top to bottom, a = 1, a = 2^ and a = 2^. The 
horizontal axis is normalized at the Nyquist frequency of the N = 64-^ distribution. 



4 DISCUSSION AND CONCLUSIONS 

We now summarize our main findings and conclusions, as 
well as indicating some directions for further study to clarify 
these issues. 

• Cosmological simulations should evidently be tested for 
discreteness effects by an appropriate, and well controlled, 
extrapolation of the relevant parameters. To recover the con- 
tinuum VP limit, we have noted that the simplest such ex- 
trapolation is to increase the particle density (i.e. ^ ^ 0) 
keeping the other relevant parameters introduced by the dis- 
creteness fixed — specifically the force smoothing e, initial 
red-shift Zi and prelC. While this may seem rather evident, 
this kind of procedure is not systematically applied in the 
literature, apart from the few isolated studies we have men- 
tioned (notably those of Melott and collaborators). More 
specifically many of the (relatively few) convergence studies 
in the literature adopt a different approach, typically de- 
creasing e in proportion to £, keeping always e ^ £. While 
such an extrapolation is not necessarily wrong, i.e., it may 
allow one to arrive at conclusions which are correct concern- 
ing discreteness effects, it has the intrinsic problem that it 
does not converge to the VP limit. Physically this means 
that such an extrapolation does not remove the non-VP ef- 
fects in the dynamics (e.g. two body coUisionality, or the 
effects described by PLT) but simply moves them to smaller 
scales. Given that the interplay of different scales in the fully 
non-linear regime of gravity is not understood, this is not a 



solid procedure. In this respect we note also that in this 
approach, additional power in the initial conditions — cor- 
responding to the extra modes which may be sampled as £ is 
decreased — is usually added. This means that structures do 
indeed form first at the smallest scales, where discreteness is 
manifestly important. Further such modification of the ini- 
tial conditions makes it difficult to identify with precision, 
as in the present study, variations which are due to dis- 
creteness. We note, however, that using wavelet techniques 
I Romeo et al. I lioOSa) have recently claimed to detect numer- 
ically discreteness effects embedded in the scatter of a set of 
cosmological simulations using different realizations of the 
initial conditions (and extrapolated power). 

• There has been some controversy in the literature about 
the widely used practice of taking results to be physical (i.e. 
representative of the VP limit) at scales below £, in simula- 
tions with e < £. We have addressed this issue with a con- 
trolled numerical study of such a simulation (with e = £/15). 
Our conclusion is that such a procedure appears to be rea- 
sonable, to a first approximation: efficient transfer of power 
from large to small scales does indeed tend to make the re- 
sults on scales below £ converge, "wiping out" the significant 
differences o n these scal es in the initial c o nditio ns (see e.g. 
' Little et al.l (|W91), and iBagla fc Prasad! l|2008i ')'l. However 
this mechanism is by no means perfect and we have demon- 
strated with our study beyond doubt that there are indeed 
measurable residual effects of discreteness at all scales, at a 
level relevant to the precision (of order a percent) now set 
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Figure 15. Normalized variance of the PS for e = 1/15, e = 1 and e = 2 at, from top to bottom, a = I, a = 2^ and a = 2"^ 



as a target for such simulations. Considering, very conserva- 
tively, only our results up to a = 2^ in Fig. [7] as indicative 
of what one would find in a typical cosmological simulation 
(i.e. starting at an initial red-shift Zi — 32), one infers a 
lower bound on discreteness effects which reach about one 
percent at the Nyquist frequency. We emphasize that these 
measures are only lower bounds, which may be very much 
below the full discreteness error. Indeed we have seen that 
in the regime (of validity of PLT) in which we can calculate 
this full error, the lower bound is (at small fc, at a = 2^) 
one order of magnitude larger than the estimated error (i.e. 
about one percent rather than the measured lower bound of 
a tenth of a percent at these scales). Most importantly the 
only way to attain greater precision, and indeed the only way 
to firmly establish the convergence to the physical limit, is to 
extrapolate to ^ ^ e. Thus, while the "common wisdom" 
is probably reasonable for the modest precision required for 
many uses of the results of these simulations, the criticisms 
formulated by some groups (notably Melott and collabo- 
rators) are fundamentally correct and further, relevant for 
the levels of precision required for some applications (e.g. 
future weak lensing observations). In this respect we note 
also that we have analysed here solely two point properties 
(essentially the PS), while Melott has emphasized that the 
numerically measured effects of discreteness are more im- 
portant in other (phase-sensitive) quantities. The methods 
used here to establish "calibrated" lower bounds on discrete- 
ness error can easily be generalized to study such quantities. 
Such a study, as well as more extended numerical studies of 



co ntrolled extrapolatio ns to the regime £ <^ e like those 
of lSplinter eralT (|l998tK using poss i bly also the methods of 
analysis employed in Romeo et al.l (|2008al ). would provide 
further insight into these issues. 



• An important element in our numerical study is the 
use of the PLT formalism. It allows us to fit analytically 
the measured dispersion in results for the PS (or, in princi- 
ple, any quantity) due to discreteness, at sufficiently early 
times. This allows us not only to "calibrate" our numeri- 
cal results, establishing that the method does really indeed 
measure discreteness eflects (rather than other numerical or 
finite-size effects) , but also gives us an understanding of the 
physical origin of these effects: a finite sampling of a fluc- 
tuation modifies its evolution with respect to the smooth 
limit. This is a physical effect of discreteness which has not 
been previously envisaged, and it illustrates very clearly that 
the widely made assumption that the effects of discreteness 
are solely those which arise from (i) missing initial power, 
and (ii) two-body coUisionality, is indeed just an assump- 
tion, which can at best be approximately correct. Indeed 
PLT describes explicitly the effect of small scales on larger 
scales which, albeit not the dominant one in the evolution 
of the gravitational clustering, is not zero when the ratio 
of these scales is finite. Such effects at large scales (i.e. sig- 
nificantly larger than I) have until now escaped detection 
in cosmological N-body simulations, even in studies which 
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looked for therrFI. Further we have noted that our results 
indicate that, apart from the very early non-linear evolution 
which "fills in" the missing power at large k, the discrete- 
ness errors at any scale continue to grow monotonically in 
time, as in PLT, throughout the whole simulation. Such be- 
haviour would naturally be explained if the physical effects 
of PLT continue to act in the non-linear regime, and indeed 
it is very plausible that this should be the case: one would 
expect that the evolution at any scale will be affected by the 
discreteness of the sampling, as in PLT, even if this sampling 
is not uniform in space as in PLT. We underline, however, 
that understanding of discreteness effects in the fully non- 
linear regime is completely lacking, and it is quite possible 
that other effects also come into play. For this reason alone 
it is important that carefully controlled extrapolations are 
systematically undertaken. 

• We have seen also in our numerical study that PLT pro- 
vides an excellent fit to the evolved power at a wavenumber 
k, until the time that this wavenumber goes non-linear, and 
indeed describes the differences between simulations on dif- 
ferent preIC for even slightly longer. This extends its valid- 
ity cons iderably beyond tha t established by the numerical 
study in iMarcos et a"l](|2006h . which showed that it extended 
only to the time when the typical relative displacement of 
nearest neighbour particles becomes of order the interparti- 
cle distance £. While this is what is expected from a naive 
analysis of the validity of PLT — requiring that the lin- 
earization in the relative displacements of the force be valid 
— it is not in fact surprising that its regime of validity ex- 
tends to the non-linearity of any given mode: to obtain a 
good approximation to the evolution of the displacement 
fields at a given scale the breakdown of PLT in describing 
the force due to particles at smaller scales is not relevant. 
The regime of validity observed is what results if one as- 
sumes that one needs the PLT linearization of the force on 
a particle to be valid only for particles at separations of or- 
der k~^ or larger. The fact that PLT does even better in 
tracing the differences between evolution from identical ini- 
tial conditions sampled on different preIC than in following 
the full evolution on an individual preIC indicates that the 
leading non-linear corrections have discreteness corrections 
which are smaller than those in PLT at linear order. A full 
study of the extension of PLT to next order (i.e. to second 
order in the Taylor expansion of the forces) should be able 
to explain this behaviour in detail. More generally, we un- 
derline that the success of PLT in fitting analytically the 
quantities we have measured shows that it can be a very 
useful instrument for controlling analytically the results of 
numerical simulations. Indeed, to our knowledge, the data 
in Fig. [To] are by far the most stringent analytic controls 
which have been placed on an A'^-body code, showing that 
GADGET can trace correctly, to a precision of as great as 
one in a thousand, differences in the PS from slightly differ- 
ent initial conditions. Thus, interestingly, the measurement 
of discreteness effects in simulations can be seen as a way 
of controlling the numerical accuracy of codes. Indeed, in 
cosmological A'^-body simulation, a reasonable goal for the 
numerical accuracy of any code is that it should measure 



See, e.g.. iLittle et all lll99lh and both the r ecent studies of 
iBagIa fc Prasad! ll2008ir and iRomeo et al.l ll2008al') . 



such effects, as it is not of physical interest to do better than 
reach this level of systematic error in the A-body method. 

• The numerical study presented was for the case of an 
initial power law PS P{k) oc fc" with exponent n — —2. We 
have also analysed fully the cases n — and n = 2, for which, 
starting from similar amplitudes of fluctuations at the scale 
£ with the same number of particles, the range of a prior 
to that at which the box goes non-linear is much greater. 
We have observed qualitatively the same behaviours, and 
in particular, the monotonic growth of the measured lower 
bounds on discreteness as a function of a. The method can 
of course be used for any initial conditions, and in partic- 
ular for the current standard ACDM model. The precise 
results for this case will depend of course, in particular, on 
what physical scale is identified with £. The use of PLT as a 
"calibrator" in this case would require its generalisation to 
this cosmology, which, as noted in Joyc e fc Marcos (2007a|) 
should be str a ightfor ward. We note that the recent study by 
iRomeo et al.l l|2008ah of this case reaches conclusions very 
consistent with those found here (and those of Melott et al 
over a decade ago) : using a wavelet analysis of a set of simu- 
lations a positive detection of discreteness errors is made for 
spatial scales smaller than of order the interparticle spac- 
ing. It would be interesting to combine in future studies 
these methods of numerical analysis with the analytical and 
numerical methods used here. 

• We have considered only numerical simulations with 
fixed £ and e, and our conclusions are valid of course there- 
fore only for this case (i .e. PM or P^M simu lations). One 
possibility, discussed bv iRomTO et all l|2008ah in their con- 
clusions, and briefly by Melott in a comment l|Melottl l|2008l ). 
see also the reply of lRomeo et al.l l|2008bl )) on this paper, is 
that the intrinsic limitations on accuracy imposed by dis- 
creteness might be addressed with numerical efficiency us- 
ing AMR type codes, with the mesh defining the resolu- 
tion of the force (i.e. effectively e) being adapted in higher 
density regions so that the condition that the number of 
particles per cell is always significantly larger than unity. 
Therefore, the idea is, one would have always a local in- 
terparticle distance smaller than the effective force resolu- 
tion scale, thus satisfying locally the condition apparently 
necessary to control discreteness effects {£ <C e) while al- 
lowing a greater spatial resolution, in denser regions, than 
that fixed by the interparticle distance £ of the initial grid. 
While such an approach would be expected to reduce greatly 
certain physical effects of discreteness — specifically any ef- 
fects due to deviations from the mean field force acting on 
particles due to particles in their immediate neighbourhood 
(e.g. by two body collisions) — our findings here lead to 
us be very cautious about this conclusion about AMR: we 
have emphasized that the discreteness effects which we have 
been able to understand physically and quantify here (us- 
ing the PLT formalism) are dynamical effects induced at any 
scale by the coupling to smaller scales at which particle sam- 
pling noise becomes dominant. When the smoothing scale is 
changed one does not undo these effects, but simply modi- 
fies them by modifying the evolution of the fluctuations at 
s mall scales. Indeed in PL T, as has been shown explicitly 
in Ijovce fc Marco3 l|2007al ). increasing the force resolution 
scale e at fixed £ does not make the evolution of the N-body 
system approximate better the physical limit. Put simply, the 
only way to reduce these kinds of effects of discreteness at 
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any given scale is to increase the particle density. Thus we 
do not consider that it is clear, in general, that an AMR type 
code can give a more accurate result (i.e. closer to the phys- 
ical model) than a standard P^M code (with e <t. £) when 
both codes use the same particle number. On the other hand, 
we would expect that an AMR code may indeed do better for 
many quantities than a simple PM code (with an effective 
e > £) at th e same particle density. In any case, as remarked 
bv iMelotd 11)08), careful tests of this or any alternative 
strategy to reduce discreteness effects should themselves of 
course be subjected to controlled tests for convergence. 
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APPENDIX A: DETAILS OF NUMERICAL 
INTEGRATIONS 



Listed in Tables lAll and IA2I are the parameter values in 
GADGET2 we have used in our "low resolution" (LR) and 
"high resolution" (HR) runs. 

GADGET2 uses adaptative time-steps, which are cho- 
sen, for each particle, using the formula: 



At = min 



2rie 



1/2 



(Al) 



where Atmax=MaxSizeTimestep, rj =ErrTolIntAccuracy 
and |a| is the acceleration of the particle in the previ- 
ous time-step (and e is the softening length). In our runs 
we have chosen the parameter MaxSizeTimestep sufficiently 
large compared to ErrTolIntAccuracy so that in practice 
only the latter parameter is relevant. GADGET2 computes 
the gravitational force using a tree-mesh technique. At large 
scales, the force is computed with a PM algorithm on a 
grid with 128'^ cells. At small scales, the force is computed 
using a tree algorithm. A node is opened (i.e. the force be- 
tween a particle and a node of the tree is computed using 
the monopole moment of the gravitational force) if 



Mf > lal/, 



(A2) 



Accuracy of time integration 


ErrTolIntAccuracy 


0.025 


MaxRMSDisplacementFac 


0.2 


CourantFac 


0.15 


MaxSizeTimestep 


0.025 


MinSizeTimestep 


0.0 


Tree algorithm and force accuracy 


ErrTolTheta 


0.7 


TypeOfOpeningCr iter ion 


1 


ErrTolForceAcc 


0.005 


TreeDomainUpdateFrequency 


0.1 


Softening lenght 


Sof teningHalo 


0.00037202380952381 


Table Al. Numerical parameters 


for our "low resolution" ri 


Accuracy of time integration 


ErrTolIntAccuracy 


0.001 


MaxRMSDisplacementFac 


0.2 


CourantFac 


0.15 


MaxSizeTimestep 


0.025 


MinSizeTimestep 


0.0 


Tree algorithm and force accuracy 


ErrTolTheta 


0.7 


TypeOfOpeningCr iter ion 


1 


ErrTolForceAcc 


0.0001 


TreeDomainUpdateFrequency 


0.1 


Softening length 


Sof teningHalo 


0.00037202380952381 



Table A2. Numerical parameters for our "high resolution" runs. 



acceleration in the last time-step and a ^ErrTolForceAcc. 
We set the option TypeOf OpeningCriterion= 1 and there- 
fore the parameter ErrTolTheta is used only in the first force 
computation, and is therefore irrelevant. 

Our LR simulations use the same range of parameters 
usually used in the literature. For example, for the VIRGO 
consortium, the parameter which controls the time accuracy 
taken as ErrTolIntAccuracy= 0.01 and the one which con- 
trols the calculation of the force as ErrTolForc eAcc= 0.005 
are considered the fiducial ones Ijenkins et al. 



are the ones which are effectively used (e.g. 
(2003); Stoehr (2006)), for a softening 



19981), and 



Stoehr et al.l 



length 



[0.0002, 0.02]£. Other works (e.g. lCrocce et all l|2006l )') divide 
their runs like us in "low-resolution" and "high-resolution" 
ones. Their "low-resolution" runs have similar resolution 
than our "low-resolution" ones (ErrTolIntAccuracy= 0.025 
and ErrTolForceAcc= 0.005), but our "high resolution" 
ones take more stringent parameters (ErrTolIntAccuracy= 
0.01 and ErrTolForceAcc= 0.002), for a value of the soften- 
ing length e G [0.02, 0.4]£. 



where M is the mass of the node of extension I at a distance 
r of the particle of which the force is computed, a the total 
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APPENDIX B: DEFINITION AND 
ESTIMATION OF TWO POINT STATISTICS 

Bl Real space 

The red uced two-point correl ation function ^(r) is defined 
(see e.g. iGabrielli at aD l|2004h ) as 

|(r) = (5(r + xM(x)), (Bl) 

where (...) means ensemble average over all the possible re- 
alizations of the system. For particle distributions ^(r) has a 
singularity at r = 0, and it is therefore convenient to divide 
it as 

^^(r) = -fo(r)+e(r). (B2) 
no 

where no is the mean number density. The quantity we give 
results for in the paper, and denote by ^(r), is a direct real 
space angle-averaged estimator of ^(r): 

1 '^■^ 
^(-) + ^- nol/(.,..)iV. g^'(-)' 

where Ni (r) is the number of particles in the spherical shell 
of radii r, r + Sr, volume V{r, dr), centred on the i*^ particle 
of a subset of Nc 4: N particles randomly chosen from the 
A'' particles of the system. 

B2 Reciprocal space 

Because we consider distributions with periodic boundary 
conditions we can write the density contrast as a Fourier 
series: 

k 

with k e {(27r/L)n|n G Z^}. The coefficients 5(k) are given 
by 

5(k) = / 5(x) exp(-ik • x)d^a:. (B5) 
Jv 

The PS is defined as 

P(k) = l(|5(k)|2), (B6) 
which we estimate with 

^W-^ E (B7) 

where N{k) is the number of vectors k' considered in the 
sum. To speed up the computation we perform a sampling 
at larger k on the vectors k'. We have checked that our 
results are robust to this choice. 



